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Operational numerical forecasting: models and products 


This issue of the Meteorological Magazine is devoted to four articles on the Meteorological Office’s 
operational numerical forecasting system. At the core of the system are the 15-level numerical weather 
prediction models, and these are the subject of the first paper. The second paper deals with the equally 
important question of the analysis of observations to provide starting conditions for the numerical 
forecasts. In the third paper attention turns to the dependent oceanographic models, for surface waves 
and storm surges, which are driven by surface winds and sea-level pressure fields from the atmospheric 
models. The further processing of results from the models to derive useful output products is the subject 
of the fourth paper. These products are required both for national purposes and for distribution to other 
meteorological services in accordance with Bracknell’s international roles. 

The articles have been prepared by staff in the Central Forecasting Branch of the Meteorological 
Office, and are intended to provide a survey from an operational viewpoint. The 15-level models and 
their data assimilation schemes were developed in the Forecasting Research Branch, and the physical 
parametrizations used in the models were formulated in the Dynamical Climatology and Synoptic 
Climatology Branches. The reader is referred to the relevant reports prepared in these Research 
Branches for full details of the various formulations and techniques. 

A further four articles, dealing with evaluation and applications of the numerical forecasts, will 
appear in the next issue of the Meteorological Magazine. 

The numerical forecasting system described in this issue, built around the 15-level models, has been 
operational since September 1982. At the time of this publication the system is well established, with a 
number of earlier problems and shortcomings remedied, and is widely recognized as being highly 
successful. Nevertheless, further important improvements are currently being researched, and 
continuing advances in performance may be expected for several years to come. 
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551.509.313 


The 15-level weather prediction model 


By A. J. Gadd 
(Meteorological Office, Bracknell) 


Summary 

The 15-level model is introduced in relation to earlier numerical weather prediction models, and the operational roles of both 
the global and the fine-mesh versions are outlined. Factors involved in the design of the model are discussed. A brief description is 
given of the main features of the formulation of the model. 


1. Introduction 


After research investigations during the 1950s and early 1960s, numerical weather prediction has been 
used operationally in the Meteorological Office since 1965. A 3-level quasigeostrophic model was used 
on the KDF9 computer until 1972, when a 10-level primitive-equation model in pressure coordinates 
was introduced on the 360/195 computer. This 10-ievel model was used in two forms: for the northern 
hemisphere north of 15°N with a 300 km grid and for a region near the British Isles with a 100 km grid. 
The 10-level model was used until 1982 and during its lifetime three different integration schemes were 
used. In 1978 important changes were made in the model’s formulation: a more accurate lower 
boundary condition was introduced and the effects of radiation were included. 

In 1982 a 15-level primitive-equation model in sigma coordinates was introduced on the Cyber 205 
computer, and this model is the subject of the present paper. 

The 15-level model is used operationally in two versions. There is a global version which is required to 
fulfil the Office’s responsibilities in Defence, as a World Area Forecast Centre for civil aviation 
(Hardman 1985), in support of marine and other commercial services world-wide, and to provide 
guidance for medium-range forecasting. The horizontal resolution used is the maximum possible for the 
original configuration of the Cyber 205 computer, with 1.5° spacing of grid points in latitude and 1.875° 
spacing in longitude. Global forecasts to 6 days ahead are computed twice daily, from 00 GMT and 
12 GMT starting conditions, and with results available by 05 GMT and 17 GMT, and, after the 
necessary further processing (Francis 1985b), are distributed to other countries in Europe and 
elsewhere. 

A second, regional version of the 15-level model has been designed principally to meet the needs of 
short-range public service forecasting for the United Kingdom (Hunt 1985) and also to support marine 
services for the European continental shelf and the Mediterranean Sea (Ephraums 1985). For these 
purposes the predictions of precipitation, surface wind and the positions of fronts are crucially 
important, and these predictions are known to benefit from a closer spacing of the grid points. This is 
achieved within the available computing capacity by using a fine mesh of points (with spacing 0.75° by 
0.9375°) over a limited area (80°W-40°E, 30°N-80°N). At the lateral boundaries of this region the 
calculations make use of information provided from integrations on the global grid. Forecasts to 36 
hours ahead are computed twice daily, with results available by 03 GMT and 15 GMT, and are 
distributed internationally for use in most of the nations within the region of coverage. The surface wind 
forecasts are particularly valuable for input to wave and surge prediction models (Francis 1985a). 

The progression of models over the years reflects the scientific advances in numerical weather 
prediction and the increasing capacity of the computers available. Thus the 15-level model is based on 
what are currently judged to be the best available treatments of the relevant phys’cs and mathematics, 
whilst its operational use is made possible by the availability of the powerful Cyber 205 computer. 
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Along with the progression of modelling techniques there has been an expansion of the area of 
coverage. This process was completed within the implementation of the global version of the 15-level 
model. However, as noted above, regional modelling is still important in that it permits a higher 
horizontal resolution, in a limited area of special interest, for a given computing capacity. 

A parallel historical development may be noted in the methods used to provide initial data for the 
numerical models. The 3-level model required objective analysis of the height field. The 10-level model 
required objective analysis of the height and humidity fields, but its initial wind values were diagnosed 
from the height analysis using the balance and omega equations. The 15-level model uses analysis and 
data assimilation of potential temperature, surface pressure, relative humidity and wind component 
data derived from observations to obtain initial fields of these variables in sigma coordinates (see the 
companion paper by Atkins and Woodage 1985). This system has its origins in work carried out as part 
of the Office’s contribution to the Global Weather Experiment in 1979. 

The increasing skill of numerical forecasts during the course of this historical development may be 
seen in sequences of verification statistics (Flood 1985). 


2. The design of the 15-level model 


- The design of the 15-level model was based on research which compared alternative techniques 
(Cullen et a/. 1981). During the lifetime of the 10-level model, three different methods for the numerical 
integration of the primitive equations had been used (Gadd 1978a), each one improving on the efficiency 
of its predecessor. Changes had also been introduced to improve the accuracy of the finite-difference 
approximations (Gadd 1978b). It was found that the numerical techniques used in the 10-level model at 


the end of its development were more efficient than, and equal in accuracy to, the alternatives and so 
were largely retained. On the other hand the representations of precipitation, convection, boundary- 
layer turbulence and radiation for the 15-level model were derived from earlier work in the Office on 
general circulation modelling and climate research (Lyne and Rowntree 1976, Richards 1980, Corby 
et al. 1977). The original formulation of the 15-level model was described in full by Dickinson and 
Temperton (1984); further improvements are constantly being researched and are introduced 
operationally as soon as they are sufficiently proven. 

At a time when several global modelling centres around the world have changed over to spectral 
models, some comment is appropriate on the retention of finite-difference methods in the 15-level 
model. The choice between spectral and finite-difference models is essentially an economic one, there 
being no intrinsic difference in the quality of result that may be obtained with either approach. In some 
centres greater accuracy has been obtained in synoptic-scale forecasting with a given computing power 
by using a spectral model. However, the efficiency and accuracy of the numerical techniques used in the 
15-level model are such that it is unlikely that any economic advantage could be obtained using a 
spectral version. Also important in the context of the Meteorological Office’s requirements is that 
exactly the same finite-difference formulation, indeed the same coding in large measure, can be used in 
both the global model and the limited area fine-mesh model. Finally, it should be noted that the balance 
of economic advantage shifts towards finite-difference models when it comes to sub-synoptic features, 
especially fronts, and this may become more important with the trend to higher and higher resolutions in 
global models. 


3. The grids 


In the vertical the sigma (o) coordinate system is used, where a is pressure (p) normalized by surface 
pressure (p,). Thus the coordinate surfaces follow the shape of the model’s orography (grid-scale mean 
values). The distribution of the 15 levels is shown in Fig. 1. In the mid-troposphere the vertical resolution 
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is about 100 mb, but near jet-stream levels there is a higher resolution of about 60 mb. A high vertical 
resolution is provided in the atmospheric boundary layer, with the lowest level approximately 25 m 
above the surface. The uppermost level is about 25 km above the surface. The 15 levels, at which all the 
model variables except vertical velocity are calculated, may be regarded as the midpoints (by log o) of 
15 layers of air. The vertical velocity (6 = Do/Dz?) is calculated at the layer boundaries. The boundary 
condition 6 = 0 is applied at the earth’s surface (o= 1) and at the upper boundary of the uppermost layer. 

The horizontal grids are defined using spherical polar coordinates ¢ (latitude) and A (longitude). With 
grid spacing Ad = 1.5° and AA = 1.875°, the global model has 121 points from pole to pole and 192 
points around each latitude circle. The fine-mesh model has a double resolution over its limited area of 
coverage. In both cases there is a staggering of model variables in the horizontal in an arrangement 
known as the B-grid and illustrated in Fig. 2. Thus the 121 X 192 (= 23 232) points of the global grid form 
the corners of grid boxes. At these corner points the surface pressure (p,), potential temperature (6) and 
humidity mixing ratio (g) are predicted, whilst the geopotential (®) and vertical velocity (a) may be 
diagnosed. The horizontal velocity components (u, v) are predicted at the centres of the grid boxes. 
Special procedures are adopted at the poles. 
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Figure 1. Sigma levels of the 15-level model. Figure 2. The arrangement of variables on the B-grid. 


4. Dynamical and numerical formulation 


The primitive equations are used in advective form, and are integrated using the split explicit 
integration scheme. This includes a modified Lax-Wendroff method for horizontal advection and a 
forward—backward scheme, with trapezoidal integration of the Coriolis terms, for the gravity-inertia 
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adjustment stage. The time-step for the integration of the advection terms is three times that for 
integration of the gravity-inertia terms. The forward-backward scheme includes special features to 
prevent the spatial lattice separation that can be a problem with the B-grid. 

In high latitudes the small east-west separation of the grid points would impose unacceptable 
restrictions on the time-steps allowed for computational stability if special measures were not taken. 
Thus stability filtering is applied as necessary in the east-west direction to allow the use of time-steps 
(15 minutes in the global model and 7//, minutes in the fine-mesh model) governed by the north-south 
grid spacings. In the global model, Fourier analysis is used around latitude circles and short-wavelength 
components of the fields are damped using coefficients calculated for each line of latitude and for each 
forecast, taking account of the maximum wind speed present. Thus the region of filtering extends further 
towards the equator in a winter hemisphere than in the corresponding summer hemisphere. Where 
stability filtering is needed in the fine-mesh model, a multipoint filter (25 points) is applied along lines of 
latitude to provide an approximation to Fourier damping. 


5. Parametrizations 


As in all modern numerical models of the atmosphere, the parametrization of physical and sub-grid- 
scale processes is an important aspect of the formulation of the 15-level model. In principle, the 
parametrizations currently available are used in the same way both in the global and in the fine-mesh 
versions of the model. However, it is sometimes necessary or desirable to introduce changes first in one 
version and later in the other. 

The current parametrizations are as follows: 

(i) Large-scale precipitation, including the ice phase and the evaporation of precipitation as it falls. 

(ii) Convection, producing a redistribution of heat and moisture in the vertical (and convective 
precipitation where appropriate) and based on parcel theory, modified by entrainment, and 
detrainment. 

(iii) Surface exchanges and boundary-layer turbulence, exchanging heat, moisture and momentum 
among the four lowest model layers and with the surface. 

(iv) Radiation, including the diurnal cycle of solar radiation, and interactive with temperature and 
humidity and with cloud inferred from humidity. 

(v) Gravity wave drag, applied at upper levels, but calculated from low-level wind and static stability 
and the variance of orographic height within each grid box. 

(vi) Horizontal eddy diffusion. 

(vii) Vertical eddy diffusion, at low latitudes only. 

The calculations of surface exchanges require values of a number of surface parameters. The 
distributions of sea-ice and snow cover are specified climatologically. Over open sea the surface contact 
temperature is analysed every day using ship, buoy and satellite measurements, and is held constant 
during each forecast. Over land and sea-ice the surface contact temperature is predicted using calculated 
components of the surface energy balance. Specifications are also required of albedo, soil moisture and 
effective thermal capacity. 


6. Lateral boundary conditions for the fine-mesh model 

The advantages available through the higher horizontal resolution of the fine-mesh version of the 
15-level model are contingent on an adequate method for updating the lateral boundaries during the 
course of the forecasts. Theoretical investigation of this problem for the primitive equations has shown 
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that a mathematically satisfactory method is very difficult to devise. Fortunately a number of pragmatic 
schemes exist which work adequately well in operational practice, and one of these is used for the 
15-level model. 

The method that is used involves a complete overspecification of the boundary conditions in that the 
tendencies of all predicted variables at all boundary points are externally prescribed. These tendencies 
are derived from values output at hourly intervals from suitable global or hemispheric integrations on 
the coarser grid and are used to update the fine-mesh boundaries every time-step. An enhanced diffusive 
damping is used at all grid points within 3 gridlengths of a boundary. Within this boundary zone the 
coefficients of the non-linear horizontal diffusion formulation are enhanced by a factor of 3 compared 
with their values in the interior regions. 
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551.501:551.509.313 
Observations and data assimilation 


By Margaret J. Atkins and Margaret J. Woodage 
(Meteorological Office, Bracknell) 


Summary 

Initial fields for the 15-level model are provided by a data assimilation scheme. Conventional and satellite observations of 
relevant meteorological quantities are processed to provide ‘observed’ values of model variables. Global analyses are produced at 
6-hourly intervals, 00, 06, 12 and 18 GMT, using data valid within 3 hours of each of these times. This is achieved by assimilating 
the processed observations into the 6-hour forecast from the previous analysis, using univariate optimum interpolation to provide 
corrections at the model grid points and additional damping to control unwanted high-frequency oscillations. Initial fields for the 
regional (fine-mesh) model are provided by a fine-mesh assimilation scheme. Analyses are produced at 3-hourly intervals by 
assimilating observations valid within 1'/, hours of the analysis time into the regional model. The fine-mesh assimilation is 
restarted every 12 hours by interpolation from a global assimilation. 


1. Introduction 


~ A necessary condition for a numerical model of the atmosphere to produce a good forecast is that the 
initial values of the model variables shold be Welt specifred. They should properly represent all the 
features, which can be resolved by the model, present in the atmosphere on a given occasion, and they 
should not contain imbalances between the variables which would excite spurious oscillations in the 
subsequent forecast. At the Meteorological Office a data assimilation scheme is used to provide initial 


conditions for both the global and regional versions of the 15-level model described in the companion 
paper by Gadd (1985). The scheme has been designed to use as much as possible of all available data, 
both conventional observations reporting at synoptic hours and satellite and aircraft observations 
which may report at any time. Before observations can be assimilated into the model they have to 
undergo quality control and be processed into a suitable form. 


2. Observations 


Meteorological data consist of measured or derived variables at different levels in the atmosphere, 
recorded at various times of the day and night. Séveral different kinds of observation are reported: 
traditional types are the surface and radiosonde reports which give excellent data coverage at synoptic 
hours over the land masses and major shipping routes. However, there are some parts of the northern 
hemisphere and vast areas of the tropics and southern hemisphere where these observations are not 
available, and it is in such regions that the satellite data have been most beneficial, providing asynoptic 
observations evenly distributed over the whole globe. Also available, along the major airline routes, are 
aircraft reports which are particularly useful to the United Kingdom in providing data across the 
Atlantic. 

All observations are transmitted in coded form on the Global Telecommunication System (GTS) and 
are received at Bracknell by the Telecommunications Branch. From there they are routed to other users 
(e.g. outstations), as well as being passed to the Meteorological Office IBM computer where they are 
decoded and stored in the Synoptic Data Bank. At this stage some basic quality control is performed 
(such as checking reported values against climatological limits and testing radiosonde ascents for 
internal consistency) and an indicator or ‘flag’ is set against any datum suspected of being in error. As 
the observations come in they are monitored by forecasters in the Central Forecasting Office (CFO) who 
examine the data together with departures from model fields using visual display units. The forecasters 
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also analyse charts of plotted data at four main levels (surface, 500 mb, 250 mb and 100 mb), and these 
are considered to be the best available estimate of the state of the atmosphere at the synoptic hours. At 
this stage forecasters are able to ‘intervene’, rejecting any data suspected of being erroneous, and 
inserting ‘bogus’ or artificial data in areas where observations are sparse and model fields are considered 
to be poor. 

The data assimilation scheme will be described in detail in the second part of this paper, but it is 
necessary to note here that it requires values of potential temperature (0), westerly and southerly wind 
components (u and v), relative humidity (UV) and model surface pressure (p, ). All observed data must be 
presented in the form of these variables in order to be used. Some types of data require only a simple 
change of variable using a standard formula: examples are the conversion of vector winds to wind 
components, the conversion of dew-point and temperature to relative humidity, and the conversion of 
temperature to potential temperature. However, other calculations could be carried out in a variety of 
ways using different approximations and assumptions, and each method gives a slightly different result. 
It is very important that the data passed to the model should be as accurate as possible, and that no 
biases should be introduced. It is equally important that the pre-processing should be carried out 
efficiently, as any time wasted at this early stage in the operational suite will cause unnecessary delay of 
the forecast products. A satisfactory balance must be found between the cost incurred and the accuracy 


attained in choosing methods of pre-processing the data. The methods used for the Meteorological 
Office 15-level model will now be described. 


Surface observations 


Surface observations comprise reports from land stations, ships and drifting buoys. These typically 
provide values of surface temperature, dew-point, wind, and mean-sea-level (msl) pressure. Exceptions 
to this are drifting buoys, which report msl pressure only, and high-level !and stations where the 
derivation of msl pressure from the (measured) station-level pressure would be unreliable. In this case 
other variables are reported, such 1s station-level pressure or the geopotential of a standard pressure 
level. Surface temperature data are not used in the model, as local anomalies are likely to render them 
unrepresentative of the atmosphere on the scale of the model grid. For the same reason, surface wind 
and humidity data are not used over land, although they are used over the sea. The most important data 
derived from surface observations are values of model-surface pressure (p,). Where the model 
orography has the same geopotential as that of the reported pressure (e.g. over the sea), no calculation 
isnecessary asp, is equal to the reported pressure. However, for most land stations some extrapolation is 
required, and assumptions must be made about the temperature structure in the layer in order to derive 
p,,- The best that can be done is to use model values, and ensure that the method adopted is consistent 
with that used to convert model p, values to msl pressure for chart output. The calculation is done using 
the hydrostatic equation integrated assuming a constant lapse rate within the layer. The lapse rate 
chosen is the ICAO (International Civil Aviation Organization) standard value (6.5 K km ') and a 
reference temperature is obtained at the model surface by extrapolating from a level above the model 
boundary layer (about 800 mb). This method ensures that where msl pressure data are equal to the 
model output values, no change will be implied in p, when the data are analysed. 


Radiosonde data 


Radiosondes provide values of temperature, dew-point and wind at significant levels of the ascent, 
together with geopotentials of the standard pressure levels. The number of levels reported varies 
according to the complexity of the atmosphere and the height attained by the sonde balloon, but is 
typically between 40 and 80. The temperature, wind and humidity data are too numerous to be passed to 
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the analysis at the reported levels, as the vertical detail could not be resolved by the model. It is therefore 
necessary to condense each ascent to a representative set of data at levels approximating to those of the 
model. This is achieved by using all reported values, weighted linearly in log pressure, to calculate mean 
values of each variable over layers centred on model levels. This has the effect of smoothing out extreme 
features such as maximum winds and tropopause temperatures which cannot be explicitly resolved by 
the model. The method is consistent with the assumptions made when output charts of standard-level 
values are produced, and ensures that geopotential (the main variable used for assessment by 
forecasters) is accurately represented. 

The reported values of standard-level geopotential are not used directly in the model, but for some 
years statistics have been collected on the differences between observed and analysed 100 mb geo- 
potentials. These are used to identify systematic differences between sondes and to derive appropriate 
corrections which are given in terms of 100 mb geopotential. In order to use this information in the 


model, equivalent temperature corrections are derived which are added to the o-level values passed to 
the analysis. 


Satellite temperature soundings 


Temperature soundings are derived from radiance data measured by the polar orbiting satellites. 
Temperature and humidity mixing ratio are the variables initially derived, but these are converted to 
thickness (or layer mean virtual temperature) and precipitable water content (PWC) for standard layers, 
before they are transmitted on the GTS. Since the Meteorological Office model does not analyse these 
variables, it is necessary to invert the calculations carried out before transmission to reconstruct the 
original data as accurately as possible. This is done by converting to layer mean temperatures, using the 
PWC data to make the correction from virtual to true temperatures. The PWC data are not used directly 
as humidity observations as the errors are rather large, particularly for higher latitudes. For soundings 
made in cloudy conditions it is not possible to obtain humidity data. The method used to derive 
thicknesses for transmission in these cases is very complicated and could not be inverted operationally. 
A simple approximate method is used instead, with a relative humidity of 80% being assumed to carry 
out the conversion of virtual to true temperature. 

Although the soundings transmitted on the GTS provide excellent coverage of temperature data, they 
are currently available only at spatial intervals of about 500 km (or 250 km in a compressed format) and 
are received at Bracknell 3 hours or more after the radiance data are recorded. They are valuable for the 
global coarse-mesh model, but are of limited use for the fine-mesh model which has a resolution of about 
75 km and runs operationally 2 hours after data time. For this reason the Satellite Meteorology Branch 
has developed a system for receiving the radiance data directly from the satellite as it passes over the 
United Kingdom, and performing the processing locally to provide temperature data directly with a 
horizontal resolution of about 80 km, typically within 1 hour of being recorded. This processing system 
is known as HERMES (High-resolution Evaluation of Radiances from Meteorological Satellites) and 
experiments are being carried out to determine the best way to use these data in the models. 


Satellite cloud-track winds 


Wind data are derived from observations of movements of clouds at various levels in the atmosphere, 
recorded by the geostationary satellites which constantly view the same area of the globe. The data are 
transmitted in the form of a wind with an associated pressure, and can be used directly in the model. 
There are certain errors inherent in these data, owing mainly to the difficulty in estimating the level of the 
cloud. However, they provide valuable information in tropical areas where other observations are 
sparse. 
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Aircraft reports 


Aircraft reports contain values of wind and temperature together with the flight level, position and 
time of observation. Observations are made every 5° or 10° of longitude along the flight path, and also at 
certain specified reporting positions. Before the data are used in the model they are checked to eliminate 
any duplicate or conflicting reports from the same aircraft. Winds and temperatures are passed to the 
analysis at the pressure level equivalent to the flight level. The conversion of flight level to pressure is 
carried out using equations based on the ICAO standard atmosphere (see World Meteorological 
Organization 1966). 


Artificial observations 


Artificial or ‘bogus’ observations are generated by forecasters in CFO to change the model fields, 
generally in areas of sparse data where they may have a better idea of the true state of the atmosphere 
from continuity, experience, or from data unavailable to the model (e.g. satellite pictures). The bogus 
observations comprise values of msl pressure, geopotential, wind and humidity at the main analysed 
levels. The geopotentials must be converted to temperatures, and the method adopted is to use model 
values to calculate the thickness change over each layer implied by the bogus data. The thickness 
changes are then converted to temperature increments, which are added to the existing model 
temperature profile to create a complete bogus ascent. The bogus data are treated by the analysis in the 
same way as any other conventional observations, but are given a higher weight to encourage the 
analysis to fit them. 


3. Data assimilation 


Many of the observations described in the preceding paragraphs may be valid at any time during a 
24-hour period. The maximum number of conventional observations are still received for 00 and 
12 GMT, this being the time at which most upper-air observations are made. For this reason forecasts 
still start from an analysis valid at those times. However, significant numbers of both surface and upper- 
air observations are received for 06 and 18 GMT and reports valid at any time are received from 
satellites, aircraft and drifting buoys. Data assimilation schemes are designed to make the best use of all 
this information, in particular, to use information valid at all times of the day to produce the best 
possible analysis valid at 00 or 12 GMT from which to start a forecast. 

At the Meteorological Office, this is achieved by producing a global analysis at 6-hourly intervals, 00, 
06, 12 and 18 GMT. All observations valid within 3 hours of each of those times are used to modify a 
6-hour forecast starting from the previous analysis. Information from observations valid at all previous 
times is thus carried forward to the current analysis in a process analogous to the use of continuity charts 
in subjective analysis. The prognostic variables of the forecast model are potential temperature (0), 
westerly and southerly wind components (wu, v), humidity mixing ratio (g) and surface pressure (p, ). 
These are carried on a latitude/longitude grid in the horizontal and 15 levels in the vertical defined in 
terms of terrain-following sigma coordinates (described in the companion paper by Gadd (1985)). The 
analysis is performed directly in the model variables (except that U is used rather than g because of 
advantages for quality control) on a grid identical to that of the forecast model, except that the number 
of points per latitude circle decreases polewards of 50° to maintain a quasi-uniform spacing. 
Observations must therefore first be processed (as described in the previous section) to provide data in 
terms of the analysis variables, but not necessarily at model sigma levels. 
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Quality control 


Univariate, three-dimensional statistical interpolation (sometimes known as optimum interpolation) 
is used in the analysis to interpolate from the observation points to the grid points. However, before this 
process is carried out, observations are subjected to quality controi making use of the same technique. 
As a first step the observations are checked against the 6-hour forecast from the previous analysis, and 
when the difference exceeds a criterion dependent on the characteristic error of both the observations 
and the forecast, the observation is flagged. In the second step, each observation (including those that 
have been flagged) is checked by comparing it with a value interpolated from the surrounding unflagged 
observations, again using statistical interpolation. If the difference exceeds a second criterion, which 
depends on the characteristic error of the observation and the interpolation error (which may be 
determined from the method of statistical interpolation), the observation is rejected. Observations 
which fail the first check may be reinstated if they pass the second check, i.e. observations which differ 
significantly from the 6-hour forecast field may be retained if they are supported by other observations. 


Statistical interpolation 


Observations which pass the quality control checks (including those performed within the Synoptic 
Data Bank and the manual monitoring of the Central Forecasting Office) are used to interpolate values 
of each of the variables at the analysis grid points. Interpolation is actually carried out in terms of the 
departure of an observation from the model field to give an increment or correction to the model field at 
each of the grid points, according to equation (1), 


ae = 3S wn + Ok a a 


observed . ° ° 
where y, is an observation of one of the variables, 6,u, v, U, p, , 


ans is the model value of the same variable interpolated to the observation point, 


Ay, is the correction or observational increment for the variable w at the grid point, 
and w is the weight given to observation w, observed. 


The summation is taken over a maximum of seven observations of the given variable surrounding the 
grid point in three dimensions. The weights w, are determined by minimizing the interpolation error over 
a large number of cases. They may then be expressed as a function of the pre-specified characteristic 
errors of the observations and model field, and the 3-dimensional spatial distribution of the 
observations around the grid point. The displacement time of the observation from the analysis may be 
taken into account by modifying the observational error. 


Assimilation 


In the most common method of applying statistical interpolation, the model field is the 6-hour 
forecast and the whole of the correction Ay, is added to this field to produce the analysis. The technique 
is usually applied in a more sophisticated form than that described above so that, for example, 
observations of wind may affect a height analysis, and vice versa, through the geostrophic relationship. 
However, even when such a sophisticated multivariate scheme is used, the simple correction of a 6-hour 
forecast by an interpolated increment will excite spurious high-frequency oscillations in a subsequent 
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forecast. These oscillations are caused by small errors in the analysis which result in imbalances between 
the mass and wind field and unrealistic values of the divergent component of the wind. They produce a 
noisy forecast and, more seriously, dissipate much of the information from the observations present in 
the analysis. The usual technique for overcoming this problem is to apply some kind of initialization. 
The grid-point values of the model variables are adjusted to satisfy certain mathematical constraints 
which remove the high-frequency oscillations while at the same time ensuring that the values remain 
close to the observations. The most popular form used today is non-linear normal mode initialization. 

At the Meteorological Office an alternative approach is used. The observations are assimilated into 
the model during tic 6-hour forecast preceding the analysis time. First the interpolation weights w, are 
calculated using values of the model forecast error appropriate to a 6-hour forecast. Then at each time- 
step of the 6-hour forecast preceding the analysis, increments are calculated at each of the grid points 
from differences between the observations and the current model state using the predetermined weights, 
by means of equation (1). However, only a small proportion (A) of these increments is added to the 
model field. This process is represented in equations (2) to (4). Let one time-step of the model be 
represented by equation (2) 


Oe. ea ck ke ae a 


where ¥ is now a vector representing all the model variables and M is an operator representing one time- 
step of the model. ¥* , ,, is the value of the forecast without assimilating any data. It is then modified by 
the observations according to equation (3), 


a ee | 


where AY, , ,, is a vector representing increments at the grid points interpolated from the observations 
for all variables at time ¢t + At. For any variable, w, at a grid point 


ica eee eR ke ee ee BE 


Equation (4) is the same as equation (1) with a slight change of notation. increases linearly from 0 to 
0.35 during the 6-hour period, so that only small changes are made to the forecast fields. This process is 
illustrated schematically in Fig. 1. For humidity, the model field (Y*, , ,,) is converted from humidity 
mixing ratio to relative humidity before the increments are calculated. After all variables have been 
incremented, the resulting field(V, , ,.) is converted back to humidity mixing ratio. Additional damping 
is included in the forecast model to control any unwanted oscillations which are generated. Since the 
interpolation is purely univariate, changes to the surface pressure and temperature fields do not directly 
affect the winds, and changes to the surface pressure field do not directly affect the temperature structure 
aloft. To assist the assimilation of the data, additional geostrophic wind increments are calculated from 
the surface pressure and potential temperature increments, and added to the wind field at each time-step; 
also, additional potential temperature increments are calculated from the surface pressure increments 
using the hydrostatic equation and added to the potential temperature field at each time-step. 

The ‘assimilation cycle’ for the global model may be summarized as follows. Observations valid 
within 3 hours of 0600 GMT are assimilated into the 6-hour forecast starting at 0000 GMT and the 
resulting field is the 0600 GMT analysis. Observations valid within 3 hours of 1200 GMT are then 
assimilated into the 6-hour forecast starting from this 0600 GMT analysis. The process continues in a 
similar way for the remainder of the 24-hour period. In practice the situation is a little more complicated 
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as forecasts have to be run as early as possible in order to be useful. At present, therefore, the 0000 and 
1200 GMT assimilations are actually performed twice, once to start the global model forecast at 0320 or 
1520 GMT and once as part of the assimilation cycle at about 1130 or 2330 GMT. This repeat or ‘update’ 
assimilation is run as late as possible to allow time for extra data to be received and for CFO to carry out 
intervention. The update assimilation is followed immediately by the 0600 or 1800 GMT assimilation. 


4) 
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iv model+data 
| 
| 
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| 
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Figure 1. Schematic diagram of data assimilation. The vertical coordinate V represents the state of the model or the atmosphere. 
The horizontal coordinate is the model time in hours. T is the validity time of the assimilation depicted and T-6 is the time of the 
previous assimilation. ¥°"*"™ is the state of the atmosphere as determined from observations at time T. ¥"*" is the 6-hour 


forecast from the previous assimilation at time 7-6. ¥"*"'* “" is the result of assimilating observations for time T into the 


6-hour forecast from 7-6. tf and 1+ Atare two successive time-steps of the forecast and ¥*,, ,,, AV, ,,, and A are as defined in 
the text. Note that ¥"*"'* ““* will not equal ¥°"""™ as allowance is made for errors in the observations. 


Fine-mesh assimilation 


Initial fields for the regional fine-mesh model are provided by a fine-mesh assimilation. This is based 
on the same principle as the global assimilation except that the assimilation is performed by the regional 
model at three-hourly intervals instead of six and, for each assimilation, observations within 1'/, hours of 
the analysis time are used. Global information from observations made at earlier times is incorporated 
by restarting every 12 hours by interpolation from a global update assimilation, and performing 12 
hours of fine-mesh assimilation at 3-hourly intervals to provide an initial field for the fine-mesh madel 
valid at 0000 or 1200 GMT. Lateral boundary conditions for this fine-mesh assimilation are provided by 
a 12-hour coarse-mesh forecast from the global update assimilation. The fine-mesh assimilation scheme 
has therefore both a higher spatial and temporal resolution than the global assimilation scheme and this 
enables more detailed analyses to be made. 
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Sea surface wave and storm surge models 


By P.E. Francis 
(Meteorological Office, Bracknell) 


Summary 

Numerical models for use as guidance in oceanographic forecasting have been used in the Meteorological Office for many years. 
The development of the present state of the art is briefly traced, with reference to the key requirements of forecasting for coastal 
defences and the North Sea oil industry. The present mode of use, mathematical formulation, and performance characteristics of 
the operational storm surge and surface wave models are separately discussed. The results of verification exercises illustrate that 
whereas forecasts of an acceptable accuracy are usually obtained there is still room for model development. Some possible avenues 
of improvement are discussed, including unification of the two continental shelf models into a single interactive system. 


1. Introduction 


The operational use of numerical models as an aid in forecasting for marine services has a long history 
in the Meteorological Office, especially if by the term ‘model’ we include such mathematical techniques 
as empirical equations based on statistical models of ocean-atmosphere relationships. Such empirical 
models were used in the Office for many years before the relatively recent advent of the more 
sophisticated models which require powerful computing resources in order to solve the basic physical 
equations. 

Following the disastrous coastal flooding of 1953 a Storm Tide Warning Service (STWS) was set up, 
working closely with the Central Forecasting Office (Synnott 1964). A system of empirical forecasting 
equations, still in use by the STWS, was first investigated for operational use during the mid 1960s 
(Keers 1966); before that time much research into the technique had been carried out at the Liverpool 
Tidal Institute (now Institute of Oceanographic Sciences) and at other oceanographic research centres. 
These equations, which relate surge residual height (the difference between actual sea level and the 
astronomical tide) to such meteorological variables as surface pressure and wind speed, were refined 
during the late 1960s (Keers 1968) and early 1970s (Hunt 1972) by the inclusion of more data and by 
investigation into methods of representing the non-linear behaviour of the surge in the shallower waters 
of the southern North Sea. More recently, in 1978, a numerical model that explicitly calculates surge 
residuals by means of solving the relevant dynamical and physical equations was introduced 
operationally. This model was also developed by the Institute of Oceanographic Sciences (Flather 1981) 
and has proved to be of great value to the STWS in their work. The model, its development and basic 
theory, are described in more detail later in this paper. 

The forecasting of surface wave conditions has been a major feature of the services offered by the 
Office since the development of the North Sea oil industry in the 1970s. Originally, forecasters at 
London Weather Centre made use of empirical growth curves (World Meteorological Organization 
1976) to gain an estimate of likely wave heights, using forecast values of surface wind speed, direction 
and fetch. Similar methods were also in use to supply the Ship Routeing Service with forecast wave data 
for the North Atlantic. A numerical model for computing surface waves in deep water was brought into 
operational use in 1976, followed by a shallow-water version in 1977. Various improvements have been 
made to these models since that time, leading to the current status of the models which will be described 
later in this paper. 

The output of the dynamical oceanographic models discussed here is used for a wide variety of 
purposes; these are discussed and illustrated in a paper to appear in the next issue (Ephraums 1985). The 
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models are run on the CDC Cyber 205 computer, making full use of the vector processing facilities. 
Results of integrations are then transferred to an IBM front-end system where post-processing is carried 
out in order to generate data in a form suitable for output purposes. 


2. The storm surge model 


(a) General discussion 


A numerical storm surge model has been in operational use at the Meteorological Office since 1978. 
The model is run as part of the operational forecasting suite, following the high-resolution limited-area 
atmospheric model, twice a day from September to May. Effectively, forecasts for up to 33 hours ahead 
are produced since the model products are not available to the forecasters of the STWS until just after 3 
hours past the initial data time. The area covered by the model is approximately 11°W to 11°E, 46°N to 
63°N, on a latitude-longitude grid with spacing of ,° x /,°. The North Sea and near continental shelf are 
well represented at this resolution, which is slightly coarser than that of the continental shelf wave 
model. The undisturbed water depth is specified at each grid point. 

Surface wind and pressure information is written up every hour during the forecast run of the fine- 
mesh NWP (numerical weather prediction) model. These data are first of all converted to wind stress and 
effective hydrostatic pressure respectively, and then linearly interpolated from the (7/,° X '/,,°) grid of the 
NWP model to that of the surge model. The ‘surface’ wind field of the NWP model is the wind from the 
lowest level of the model, located at approximately 25 m above sea level; the implications of using a wind 
at this level are examined in a later section. 

Two runs of the surge model are performed in order to obtain the residual height data. In the first run 
the tidal predictions for each grid point are calculated using 6 tidal components. These data enable a 
second run to be made which includes the meteorological effects in addition to the tidal terms. The surge 
residuals are then obtained by subtracting the results of the purely tidal run from that which includes 
both tidal and meteorological terms. In this way the tide-surge interaction, which is of greatest 
importance in shallow water, is included in the calculations. 

The initial conditions for a surge model forecast are at present taken from the results of the previous 
12-hour forecast. In this way any errors due to a poor wind and pressure forecast are carried over and 
perpetuated. A method of avoiding such a situation is to introduce a ‘hindcast’ run, similar to that 
performed in the surface wave models (see later). Such a system has been tried in the past and proved to 
be worth while, and there are plans to introduce it for the present storm surge model. 


(b) Mathematical formulation 


The mathematical basis of the model is formed by the depth-integrated equations of motion and 
continuity. Following the most recent and comprehensive description of the operational model (Procter 
and Flather 1983), the complete equations are approximated by: 
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where the notation is as follows: 

f=2Msind _ the ‘Coriolis’ parameter; 

t, A, @ time, east-longitude and latitude respectively; 

U, v, components of the depth mean current g; 

l elevation of the sea surface; 

h undisturbed water depth; 

D=h+C total water depth; 

R,,g, geophysical constants (earth radius, earth rotation rate, gravity, density of sea water); 
P. atmospheric pressure at sea level; 

F., Gs components of wind stress 7, on the sea surface; 

Fy, Go components of stress on the sea bed, 7». 

The stress on the sea bed is calculated by assuming a quadratic stress law in order to relate the stress to 
the depth mean current, i.e. 


Tr = -0.0025 p qiq|. 


Similarly the surface stress is computed from the atmospheric wind data by means of the relationship 


Ts = pa 10-° (0.63 + 0.066 |W |)W |W | 


where p, = the density of air (1.25 kg m~*) and W is the surface wind. 

Horizontal boundary conditions are given by assuming zero flow normal to coastlines, and by seeking 
to prevent artificial reflections from the open boundaries of disturbances generated within the model. 
The latter effect is achieved by allowing outward propagation as free progressive waves, while ensuring 
the generation of the tidal constituents at the boundaries. The details of the finite difference schemes 
adopted for the numerical integrations of equations (1)-(3) may be found in Procter and Flather (1983). 
By excluding very deep water areas off the continental shelf, the time-step of the integration is allowed to 
be fixed at 144 seconds. 


(c) Performance characteristics 


During the winter season of 1983/84 the STWS found that in general the model overpredicted positive 
surge residuals, especially in the southern areas of the North Sea. The data used for verification were 
taken from reliable tide gauges located at ports around the British coast. Some error statistics for four 
locations are given in Table I. The tendency for overprediction is shown by the positive mean errors, and 
also by the larger extreme values of over-forecasts, much larger numerically than the extreme under- 
forecasts. A small improvement in forecast quality with reduction in forecast period can be seen, but this 
is really very marginal. A possible explanation of the overestimation seen in these recent verification 
figures may be found in the direct use of winds taken from the lowest level of the NWP model. 
Traditionally, winds at 10 m above sea level have been used throughout the development of the surge 
model, especially in the derivation of the surface wind stress relationship (Smith and Banke 1975). It is 
thus possible that further modelling of the wind at the output stage, using the thermal structure of the 
NWP model, in order to generate winds nominally at 10 m above sea level, would be a beneficial 
development. 


3. The surface wave model 


(a) General discussion 


The surface wave model at present run on an operational basis in the Meteorological Office is the 
result of an evolving series of such models, in use since 1976. Similarly to the surge model, the wave 
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Table I. Error statistics of storm surge model residual estimates, winter 1983/84 


Port, and time of Number Root-mean-square Mean Largest errors 
issue of forecast of cases error error overf'cast/underf’cast 
before high water metres 


Tyne -12 40 0.20 +0.09 +0.48/-0.19 
-4 33 0.18 +0.04 +0.46/-0.27 
Immingham -12 39 0.23 +0.17 +0.55/-0.11 
-6 33 0.23 +0.13 +0.55/-0.33 
Lowestoft -12 30 0.44 +0.34 +0.89/-0.24 
-6 30 0.33 +0.21 +0.90/-0.33 
Sheerness -12 36 0.38 +0.22 +1.00/-0.33 
a7 24 0.31 +0.22 +0.88/-0.11 


models (three versions) are driven by wind fields extracted from the operational NWP models. All 
models are run twice daily, but throughout the whole year rather than just in the winter season. The three 
versions in operational use cover different areas on different scales, and are run for different forecast 
periods (see Table II). The only difference in formulation between these three models is the inclusion of 
extra shallow-water terms in the continental shelf model. The model integrations are at present carried 
out on a polar stereographic grid, following the appropriate interpolation of the wind fields from the 
latitude-longitude grids of the NWP models. The wind input for the wave models is assumed in the basic 
theory to be at a height of 19.5 m, hence the use of direct lower-level winds from the NWP models is not 
inappropriate here as it probably is for the surge model. 


Table II. Physical specification of operational wave models in the Meteorological Office 


Model Resolution Origin and frequency Forecast Depth-dependent 
at 60°N of wind field data period terms 


Northern hemisphere 150 km Global NWP model, 3 days No 
every 1'/, hours 

Continental shelf 25 km Fine-mesh NWP model, 36 hours Yes 
every /, hour 

Mediterranean 50 km Fine-mesh NWP model, 36 hours No 
every , hour 


Unlike an NWP model, where there are sufficient observed data to allow an ‘analysis’ (or starting 
field) to be calculated for all the different variables, the wave model has to be initialized by means of a 
‘hindcast’. This is essentially the repeat of the first 12 hours of the previous forecast (run 12 hours ago) 
but this time using improved wind fields, improved by making use of the intervening analysis runs of the 
NWP model. There are insufficient measurements of the directional wave spectrum to allow a wave 
analysis and, additionally, the theoretical mathematical work on which to base such an analysis has yet 
to be carried out. The hindcast technique yields a series of fields, 12 hours apart, which are the best 
description of the actual wave fields that can at present be derived. The archive of these fields has many 
uses; these are discussed in a companion paper (Ephraums 1985). 


(b) Mathematical formulation 

The theoretical basis of the wave model is the representation of the wave field as a spectrum — 
essentially the variance of the sea surface elevations that are present in a statistically random fashion. 
The behaviour of individual components of the spectrum can be calculated from known theories, and 
the more commonly recognized quantities such as wave height and period can be calculated from 
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moments of the spectrum. A simple example of a directional wave spectrum is shown in Fig. 1. The 
model at present has 16 equally spaced directions and 14 irregularly distributed frequency values, 
ranging from 0.04 Hz (25 seconds period) to 0.324 Hz (~3.1 seconds period). Thus at each model grid 
point 224 discrete spectral energy density values are stored. 
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Figure 1. The limiting energy spectrum for a wind of 20m s“' from 360°. 


The variation of these values in time and space is controlled by means of integrating the energy 
balance equation: 
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where the following notation is used: 

E(f, 6, x, t) is the energy density in a spectral component with frequency /, travelling in a direction 0, ata 
location x at time #; 

Cf, 0) is the group velocity of the component (a function of depth); 

Sin represents the energy input from the atmosphere; 

Sus represents energy loss; and 

Sn represents the internal redistribution of energy due to non-linear interactions. 
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Full details of the mathematical and physical theories used in the model may be found in a 
comprehensive paper (Golding 1983). It is sufficient here to outline the processes in a more general way. 

The advection, or propagation, term is calculated by means of standard finite difference techniques, 
with appropriate treatment at coastal and open boundaries in order to prevent reflection. The refraction 
term is solved by adopting a continuous form of Snell’s Law, more familiar in optics, in order to 
calculate the effects of irregular bottom topography on the component waves. The source and sink terms 
Sin, Sus contain a number of processes which attempt to represent the currently held theories concerning 
wave growth and decay. In simple terms the growth of sea waves is due to the input of atmospheric 
energy through boundary-layer transfer processes which are still not completely understood. What is 
clearly apparent is that growth is inhibited by a limit on wave size and shape, beyond which breaking 
takes place with an accompanying loss of energy. The growth and decay processes in the model (for 
deep-water states) are balanced so that the limiting form is attained but not exceeded, and also so that 
observed growth rates are reproduced. In shallow water an extra dissipation term is introduced to 
represent the energy loss through frictional effects on the sea bed. 

Finally, the non-linear interaction term, Sn, is introduced in order that the eventual spectral shape is 
close to that observed in physical reality. It is known from an extensive series of measurements that the 
spectrum of a wind sea (i.e. that portion of the full spectral range within the influence of the local 
generating wind) is maintained in a recognized shape by means of internal energy transfers from high to 
low frequency. The calculation of the exact transfer terms is, however, at present beyond the capacity of 
most computing systems, especially if real-time forecasts are required. Hence we have adopted an 
empirical approach which assumes a spectral form, the JONSWAP spectrum (Hasselmann et a/. 1973), 
and redistributes the wind sea energy, without loss, into that form. The main difficulty with this method 
is the problem of separating wind sea from swell (i.e. that portion of the spectrum not under the local 
influence of the wind) in order to generate the spectral shape, but the model results as discussed in the 
next section confirm that the present empirical scheme is giving an acceptable approximation to a very 
complex physical process. 


(c) Performance characteristics 


The wave model can be assessed in two distinct ways: firstly by comparison of values of significant 
wave height (approximately that height assessed by visual estimation, more specifically the mean height 
of the highest one-third of the waves present) as calculated by the model and as observed and reported 
from maritime stations; and secondly by means of comparison with measured wave data, for example 
that from a wave rider buoy or similar device. The first verification technique is performed operationally 
at Bracknell, where incoming wave data from synoptic reports are routinely compared with model 
results. Table III shows some error statistics from three representative stations. Ocean Weather Ship ‘L’ 
is in deep water, the oil platform Statfjord is near the edge of the continental shelf, depth about 130 m, 
and the platform K-13 is in the much shallower regions of the southern North Sea, depth about 30 m. 
Points to note from the figures in the table are: 

(i) Errors grow with observed wave height, but not so markedly with forecast period. 

(ii) The effects of the extra friction terms in the continental shelf model are apparent in the results for 
Statfjord where, although the water depth is still appreciable, the coarser-mesh model is yielding higher 
wave height forecasts throughout all ranges of observed height. 

Details of a comparison against data from a wave rider buoy can be found in Fig. 2. The fine-mesh 
model has given very acceptable forecasts for a location about 500 m offshore in 13 m of water. Such a 


performance confirms the potential that the model has for many shallow-water applications (Ephraums 
1985). 
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Table III. Error statistics of significant wave height in operational wave models in the Meteorological 
Office, 1983 and 1984 


0-3 m 3-6m 


Station Forecast a 
time model N o N Xx 


OWS ‘L’ 00 NH 466 (0. 0.9 633 0.3 
57°N, 20°W 24 467 0. 0.9 647 0.3 
0. 
0. 


Statfjord 00 NH 586 3 0.7 375 (0.9 : 41 
61.2°N, 1.8°E 24 598 3 08 380 =6—(0.9 ‘ 41 
00 cS 580 0.1 0.6 371 = 0.6 . 41 

24 587 0.2 0.7 374 = (0.7 : 40 


K-13 00 cS 1037 0.0 0.3 86 -0.3 : 0 
53.2N, 3.2°E 24 1055 0.2 0.5 87 -0.1 ‘ ‘0 


Wave height ranges are for observations; N= number of cases; X= mean error in metres; o = standard deviation of error in metres; 
NH = northern hemisphere; CS = continental shelf. 
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Figure 2. Comparison of wave heights, actual and forecast, for Seaford, January 1984 (by courtesy of Southern Water, and 
Hydraulics Research, Wallingford). 
—________ observed wave heights from wave rider buoy 
a 7 @ 12-hour forecast (for 12 GMT made at 00 GMT) 
oO fe) © 24-hour forecast (made at 00 GMT the previous night) 
Forecasts at location 50.6°N, 0.8°E; Seaford wave rider buoy at 50.76°N, 0.08°E. 


4. Future developments 


The present situation, in which the storm surge model and the fine-mesh model are run separately, is 
far from satisfactory since a unified model would add far more to the accuracy of the results than just an 
increase in operational efficiency. The potential advantages in a unified model are threefold. Firstly, 
long-wave propagation is strongly influenced by water depth, hence tidal and surge influences on water 
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depths in areas of shallow mean water level are likely to be significant for wave forecasts. The southern 
North Sea is a particularly interesting area in this respect. Secondly, it is well known that currents have a 
marked effect on waves, and it is quite likely that the depth-integrated currents of the surge model could 
be used to give better wave forecasts in regions of strong tidal currents. Thirdly, the surface stress which 
drives the surge model is directly linked to the roughness of the sea surface. Thus a knowledge of the state 
of sea should result in more accurate computations of the stress input for surge modelling. 

Both kinds of operational oceanographic model will benefit from the assimilation of measured data 
on a routine basis. For the surge model the requirement is a balanced set of current and height data; for 
the wave models a fairly comprehensive description of the directional spectrum is necessary. Such data 
are at present not readily available and it may be some years before data in large enough amounts, e.g. 
remotely sensed from satellites, make such an assimilation process viable. In the meantime some 
research work is necessary in order to solve the mathematical difficulties inherent in the process. 

Research work into wave and surge modelling is at present coming to the fore in many European 
countries, and a joint initiative for co-ordinating further work has recently been inaugurated. The WAM 
group (WAve Modelling), with members from nearly all the maritime European nations, intends to 
address itself to the many problems that remain to be solved in order that real-time oceanographic 
forecasting for a variety of marine applications will achieve even better results than at present. 

As part of this program the Meteorological Office is currently performing experiments with a wave 
model on a latitude-longitude grid, with the intention of providing a full global forecasting capability in 
1986. At the same time the new model is to be used as a ‘test bed’ facility with which to investigate the 
progress of some of the work in the WAM initiative. 
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Output products of the Bracknell numerical weather prediction models 


By P. E. Francis 
(Meteorological Office, Bracknell) 


Summary 

The data produced by the numerical weather prediction models are further processed in order that a comprehensive set of 
output products can be generated operationally. The first step is the generation of an output data set which contains both primary 
and derived variables, arranged on fields at constant standard pressure levels. The derived variables are those which are not 
explicitly contained in the forecast model data. 

Charts and coded bulletins are then constructed from the output data set, using a variety of formats requested by users, and 
internationally agreed. Lists are provided of the current charts and code bulletins available. A discussion of back-up arrangements 
concludes the presentation of output products. 


1. Introduction 


The numerical weather prediction (NWP) models described in a companion paper (Gadd 1985) yield 
large amounts of data when run operationally. These data (for different variables, model grid points, 
model levels, forecast times) are internally arranged and represented in a form that is convenient for the 
rapid integration of a numerical model. Such a representation is, however, far from convenient for the 
production of the large quantity of output, in many forms, that is needed to satisfy the requirements of 
modern meteorological services. In order that a versatile range of output products can be easily 
generated it is first necessary to rearrange the forecast model data fields into a more suitable format. At 
the same time it is convenient to generate (by further modelling techniques) those fields which are not 
explicitly carried in the models, but which are required for operational purposes. These ‘derived’ fields 
are listed in Table I, together with the other standard formats for explicitly held variables. 

The output variables, held in a comprehensive data set, are then transferred from the Cyber 205 fast 
vector processor to a conventional ‘front-end’ processor, currently an IBM 3081. After reformatting, the 
output data set is available for access by authorized users, as well as being the source of all the 
operational output data streams. A versatile suite of access routines is available which allow data to be 


Table I. Variables stored in output formats, both basic and derived 


Model variable Basic output forms Derived output forms 


Wind Speed and direction on Vertical velocity at standard levels. 
standard pressure levels Maximum wind, height, speed, direction. 
Surface (10 m) wind, speed, direction. 
Vorticity at standard levels. Clear-air 
turbulence index (200, 250, 300 mb). 


Temperature Temperatures and heights of Tropopause temperature, pressure, 
Surface pressure standard pressure levels. height. Height and pressure of freezing 
Mean-sea-level pressure level. Surface (1.5 m) temperature. 
Wet-bulb potential temperature 
(500, 850 mb). 


Humidity Relative humidity at 
standard pressure levels 
2500 mb. 


Precipitation Dynamic and convective 
rates and accumulations. 
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extracted and then interpolated on to a wide range of latitude-longitude grids, and also on to different 
map projections, including Polar Stereographic and Mercator. Another advantage of retaining output 
products on a front-end processor is the ease of communication with other computing systems, and the 
fact that incoming back-up data from the National Meteorological Center (NMC) in Washington DC, 
and the medium-range data from the European Centre for Medium Range Weather Forecasts 
(ECMWFE) can be collected and stored in the same output format, ready for presentation in chart form 
using standard computer software. 


2. Output modelling techniques 


The first stage in the process of generating output products (Fig. 1) is the construction of a 
comprehensive data set which consists of the required basic and derived variables, stored on pressure 
levels. Table I briefly sets out the current contents of the output data set. Not all of the derived fields are 
generated by both versions of the NWP model (i.e. global and fine-mesh versions), but all of the basic 
output forms are common to both. Variables directly associated with aviation forecasting (tropopause 
and maximum wind, for instance) are generated only from global model data. 
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d-—levels 
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Front-end processor ; > : 
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Figure 1. Organization of data set for generating output products. 


(a) Wind modelling 


The initial step in the wind modelling is the calculation of the vertical velocity. This calculation is 
performed while the horizontal wind components are still on o-levels, and is carried out in a manner 
similar to that employed in the forecast model. Divergences at every grid point are calculated for each 
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model layer and summed from the top downwards. In the same manner, absolute vorticity is also 
calculated from the horizontal wind components, before interpolation on to pressure levels is carried 
out. The required interpolation, for both vertical velocity and absolute vorticity, is calculated by 
assuming that the variables change linearly with the logarithm of pressure, written log (pressure) below. 

The next step is to calculate the speed, direction and pressure level of the maximum wind. The wind 
components are first combined to give a speed and direction, then a search is made for the o-level that 
contains the maximum value of wind speed. The maximum of the implied continuous profile may then 
be assumed to lie within one model layer on either side of the designated o-level. A cubic spline function 
is fitted to the wind speed values outside of the boundary layer, and by means of straightforward but 
detailed calculus and algebra it is possible to calculate both the profile maximum value and its level. An 
enhancement of the maximum value is not at present thought to be necessary since the improved 
definition of jet streams in the present model gives acceptable results in most situations. The direction of 
the maximum wind is found by means of interpolation in log (pressure) once the level has been 
determined. 

Consistent with the vertical interpolation used for other wind parameters, wind speed and direction 
on standard levels are determined by assuming a linear variation with log (pressure) of the u, v 
components between the o-levels. The maximum wind information is used as an extra level in order to 
give consistency in that region. Standard levels that fall below model orography are assigned the wind 
speed and direction of the lowest o-level. 

Another derived ‘wind’ parameter is the clear-air turbulence (CAT) index, which is calculated using 
the formulation established for the previous NWP model (Dutton 1980). Work on updating the 
formulation, using data from the new NWP model related to aircraft reports of CAT, is now under way, 
and should lead to improvements in the definition of this very important aviation index. 

Finally in this section on wind modelling, we turn to the problem of calculating a surface wind defined 
as the wind at 10 m above the ground or sea surface. For many forecasting applications such a wind 
parameter is an essential element and it is important that realistic values are derived. The technique used 
during the output modelling calculations is consistent with that employed in the forecast model for 
representing boundary-layer effects. An integral part of the system is the calculation of a boundary-layer 
Richardson number, taken to be a function of wind, temperature and humidity values at the surface and 
lowest o-levels. This Richardson number then determines the profiles used to evaluate both surface wind 
and temperature over land and over sea. The profiles vary with the implied stability in the boundary 
layer, and contain the implicit contributions of surface roughness length and drag coefficient. The 
direction of the surface wind is taken to be that of the wind in the lowest o-level. 


(b) Temperature modelling 


The first stage in the process of calculating the heights and temperatures of standard levels is to 
convert the potential temperatures on o-levels, used in the forecast model, into the associated 
temperature values. Then the thickness between o-levels can be calculated using the o-level virtual 
temperatures. In this way the effects of humidity on the thicknesses are included. The heights of the 
o-levels can then be derived by summation of thicknesses from the lowest level upward. The temperature 
at each standard pressure level is calculated firstly by locating it between two o-levels, and then by 
assuming a constant temperature lapse rate between those levels. The height at each standard pressure 
level is calculated by assuming a constant virtual temperature for the layer between the pressure level 
and the lower (in height) of the two reference o-levels. 

For pressure levels below the model orography it is necessary to adopt a pragmatic approach in order 
to obtain acceptable temperature and height values. The associated calculation of pressure at mean sea 
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level is performed in a consistent manner. Height values below orography are primarily required so that 
realistic contour charts can be produced for the forecaster. Such charts require to be smooth, and to 
appear as if the intruding orographic feature was not present. Accordingly it is acceptable to extrapolate 
from a ‘free atmosphere’ level in order to obtain the temperature for a subterranean standard level, thus 
avoiding the effects of the boundary layer. The extrapolation uses a standard lapse rate of 0.0065 K m='. 
The height at the subterranean standard level is obtained by subtracting the thickness of the layer 
between the lowest o-level and the standard pressure level concerned, again calculated using the 
assumption of a constant virtual temperature throughout the layer. The pressure at mean sea level is 
given by an equivalent calculation using the orographic height and an extrapolated temperature profile 
as above. For zero orography the pressure at mean sea level is of course equal to the surface pressure. 

The next task is to locate the tropopause and assign a consistent temperature and height. The essence 
of the method employed to do this is the assumption that the tropopause should be associated with a 
layer mean lapse rate of less than 0.002 K nr, at a height above that of surface inversions and other low- 
level temperature features. Layer mean lapse rates are therefore calculated, from the o-level data, anda 
search is carried out (above 500 mb) for the relevant critical lapse rate. Having found the first layer for 
which the mean lapse rate is less than the critical value it is a straightforward algebraic task to calculate 
the exact o-level for which the lapse rate equals 0.002 K m-', assuming a linear variation in the lapse rate 
between the mid-points of adjacent layers. Following further algebraic manipulation the temperature of 
the tropopause can be determined, and then in the same manner as described earlier the height can be 
calculated. 

Three other temperature-dependent variables are also derived. The freezing level is located by 
searching upward through the o-level data for the first value of temperature less than 273 K. The freezing 
level height is then calculated by assuming a constant lapse rate through the appropriate model layer. 
The technique for modelling surface parameters was outlined earlier in the discussion of wind variables. 
By this method a temperature at 1.5 m is modelled and made available for output. Finally, values of wet- 
bulb potential temperature are calculated by an iterative technique after the necessary temperature and 
humidity information has been transferred on to pressure levels. 


(c) Humidity and rainfall 

The forecast model produces values of humidity mixing ratio on o-levels whereas for output purposes 
relative humidity on pressure levels is required. The first step in the transformation process is the 
calculation of saturation mixing ratio, performed by using a look-up table for saturation vapour 
pressure with reference to the o-level temperature. The look-up table assumes saturation with respect to 
water for temperatures greater than -5 °C, and with respect to ice for temperatures less than-10°C. For 
intermediate values the table has an interpolated entry between the two saturation states. The 
information is then at hand to be able to first calculate saturation mixing ratio and then relative 
humidity. Interpolation to pressure surfaces is again achieved by assuming a linear variation in 
log (pressure). 

Rainfall information in the models is available as the result of both convective and dynamic processes, 
the basic parameter being an instantaneous rate. By accumulating the rates at each time step and 
multiplying by the appropriate constant it is also possible to derive values for dynamic and convective 
accumulations over a period. This calculation is performed in the forecast model. Enhancement of the 
grid-point convective values in order to represent possible convective activity in a fraction of a grid 
square is at present performed empirically; but the use of relevant information from the convective 
parametrization in the forecast model may lead to a more satisfactory procedure for such enhancement. 
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3. Chart output 


The conventional presentation of most forecast products is by means of a chart, either containing 
isopleths of the required field or having numerical values printed on a regular grid. Chart products form 
a very important part of the output stream from the numerical models run at Bracknell. Including all 
kinds of charts, and those required for back-up purposes, nearly 1000 are produced in a 24-hour cycle. 
Not all of these are transmitted by facsimile broadcasts of course. Nearly half of them are for back-up 
purposes, i.e. for use if a later forecast run fails. Many are for use in the Central Forecasting Office and 
are not designed for wider dissemination. About 100 charts are prepared for facsimile broadcast in a 
24-hour period, including a substantial number of hand-drawn products. These latter charts are 
prepared by forecasters at Bracknell, using forecast model products as a guide, along with other 
information such as satellite pictures and later observations. The final chart represents a ‘man-machine 
mix’, that is, a synthesis of computer speed and efficiency with human experience and intuition. Both 
hard-copy and visual display unit presentations are available to the forecaster at Bracknell (Fig. 1) for 
the preparation of hand-drawn products. 

Three chart examples are shown in Fig. 2, and other examples are to be found in companion papers. 
The computer-drawn charts are produced on a Calcomp 1581 device, using a software package which 
extracts data from the comprehensive output data set (described earlier), and prepares it for either line 
drawing or grid plotting on a versatile range of chart areas and projections. Tables II and III summarize 
the charts currently available on facsimile broadcasts. The more acceptable modern method of 
disseminating forecast data is by means of coded broadcasts using higher speeds, which enable a wider 
range of data to be transmitted. This system will be discussed in the next section. 


4. Coded output 


(a) Bracknell’s dual role 


The Meteorological Office at Bracknell has a dual function as far as the dissemination of coded 
forecast data is concerned. As a Regional Meteorological Centre (RMC) for Region VI of the World 
Meteorological Organization (WMO), Bracknell makes available numerical analysis and forecast data 
covering Europe and the North Atlantic to other national meteorological centres. Acting as a World 
Area Forecast Centre (WAFC) of the World Area Forecast System (WAFS), initiated by the 
International Civil Aviation Organization (ICAO), Bracknell provides data with a global coverage, 
both directly to airlines and also to regional centres of the WAFS. The provision of these data requires 
the use of agreed code forms, so that rapid transmissions from centre to centre are achieved without 
confusion or ambiguity of meaning. 


(b) The code formats 


The rapid progress in telecommunication techniques, leading to high-speed lines and computer-to- 
computer connections, has resulted in a parallel development in the use of codes for transmitting 
numerical forecast data from one centre to another (Fig. 1). The most widely used code form is WMO 
GRID code (FM 47-V) which was devised for use on low-speed lines, using a character-based format, 
with relatively high overheads associated with possible output on line-printer devices. The majority of 
coded bulletins from Bracknell are in this format, and will remain so for the forseeable future; however, 
more efficient forms have been devised (see below) which should prove attractive enough for a major 
change in emphasis over the next five years. GRID code forms the mainstay of the network between 
national weather centres in Region VI. 


The data required for flight planning purposes by airlines have traditionally been supplied in Aviation 
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Figure 2. Examples of forecast charts. 
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Table Il. List of available chart products using Analogue Facsimile transmissions 


In the following tables all the computer-produced charts are based on data from the global numerical weather prediction model. 
Most of the subjectively drawn charts (M) are also largely based on the products of the same model or its derivatives. 


(a) Analysis products 


Reference time (GMT) Method of production 
and area of coverage fe oe 

(see Table IIT) meme _ 
Field identifier 06 12 18 


Surface pressure, with F F F M 
front and centres M 


850, 700, 500, 300, 200 y A 
and 100 mb contours, and 
500-1000 mb thickness 


Cc 


500 mb contours 


500 mb contours and 
500-1000 mb thickness 


Sea surface temperature 
Sea ice 


State of sea 


(b) Forecast products 


Validity time Method of production 
(hours after data time) eins conti an C 
and area of coverage Grid path are 
(see Table Til) Manual input —M 
24 48 72 96 120 
Surface pressure, < a 
with fronts and 


centres 8 c Cc a 


Soars 


850, 700, 500, 300, 
200 and 100 mb 
contours 


500-1000 mb thickness 


Ss 
on ZEEE 


— 
N 


ns 


Winds and temperatures 
(DDFFFTT) at 850, 700, 
500, 400, 300, 250 and 
200 mb 


=) 
aa aN 


N 


Winds and temperatures 
(DDFFFTT) at 500, 300, 
250 and 200 mb 


Significant weather 
(jets, tropopause, 
CAT, etc.) 


State of sea (wave 
height, direction) 
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Table Il. Bracknell (Regional Meteorological Centre) — list of chart areas. The projection used for all of 
these charts is Polar Stereographic. 


Code 
letter 


(48°N, 
(28°N, 
(42°N, 
(29°N, 
(S7°N, 
(69°N, 
(42°N, 
(80°N, 
(53°N, 


145°W) 
145°W) 
90°W) 
155°W) 
96°W) 
111°W) 
112°W) 
0s°w) 
50°W) 


Area 


(32°N, 68°E) 
(24°N, 54°E) 
(66°N, 90°E) 
(28°N, 63°E) 
(719N, 71°E) 
(37°N, 50°E) 
(60°N, 32°E) 
(44°N, 32°E) 
(53°N, 65°E) 


(18°N, 12°E) 
(09°N, 05°W) 
(30°N, 20°E) 
(08°N, 06°W) 
(46°N, 14°E) 
(19°N, 10°E) 
(28°N, 10°W) 
(28°N, 05°W) 
(26°N, 34°E) 


(28°N, 71°W) 
(02°N, 84°W) 
(20°N, 40°W) 
(08°N, 85°W) 
(38°N, 48°W) 
(34°N, 55°W) 
(21°N, 74°W) 
(43°N, 41°W) 
(26°N, 18°W) 


Scale in 
millionths 


1:20 


(31°N, 20°W) 
(67°N, 37°W) 


(23°N, 150°E) 
(70°N, 18°E) 


(06°S, 102°E) 
(47°N, 07°E) 


(05°S, 33°E) 
(45°N, 19°W) 


ZrALTONnMoINw>y 
SRSBBBSSSR 


Digital Forecast (ADF) code. ADF was designed for a specific purpose and thus differs from GRID 
code (and its successors) in that the data are arranged in columns rather than in horizontal fields. The 
code uses a ‘packed-decimal’ format, i.e. more efficient than GRID code, but is still character based and 
thus carries a large overhead penalty. Bracknell currently provides forecast data in this code form to 
British Airways, Scandinavian Airlines System, and the Société Internationale de Télécommunications 
Aéronautiques. 


In order to exchange the very large amounts of data required for the operation of the WAFS it was 


necessary for a much more efficient code form to be devised. As an interim measure Bracknell and the 
National Meteorological Center in Washington developed a bit-orientated code form, H code, which 
was implemented between the two centres in 1984. Subsequently WMO extended and developed the 
concept; a more widely agreed bit-orientated code, GRIB (GRIdded Binary), was ready for use at 
Bracknell from 1 January 1985, and the complete switch-over from H code took place in April 1985. 
This bit-orientated code has been designed for use on fully automated computer-to-computer systems 
and is compatible with the more efficient telecommunication protocols required for transmission on 
high-speed lines. The supply of data from Bracknell to automated systems (OASYS) at London Weather 
Centre and Heathrow Airport is an example of this computer-to-computer facility. The major part of 
the Bracknell coded products list is now available in GRIB code. 


(c) Products available in code 


Table IV contains a concise summary of the products that are currently available in coded form from 
Bracknell. Three distinct classes of product are identified: those taken from the fine-mesh NWP model, 
primarily for use in WMO Region VI; regional products from the global coarse-mesh NWP model, again 
for use in Region VI and extending into the medium range; and global products from the coarse-mesh 
model, formulated primarily for use in the WAFS but also available for supply to other interested users. 
Details of the elements contained in the bulletins, areas covered, bulletin resolution and forecast times 
are all to be found in Table IV. 

The most recent additions to this list of products are the very-high-resolution (1¥,° X 1,°, every 3 
hours) bulletins taken from fine-mesh NWP model data. The Royal Netherlands Meteorological 
Institute requested such data for input to a scheme being developed for forecasting of maritime weather 
and state of sea conditions along the Dutch coast. Other recent requests for data, which have been met 
by bulletins from the extensive range available at Bracknell, came from Zimbabwe, Mauritius, Hong 
Kong and the Seychelles. 
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Table IV. List of available forecast fields from Bracknell. 


Forecast runs are made from both 00 and 12 GMT data time. Bulletins are coded in WMO GRID (FM 47-V). Details of bulletin 


headings and contents can be supplied on request, including details of GRIB-code bulletins. 


Elements 


Resolution 


Forecast 
times 


Regional fine-mesh products 


Surface pressure, surface 
wind (10 m), 
accumulated 
precipitation. 


Heights, winds and 
temperatures for 850 
and 500 mb. 


The 3 areas all span 
32Y,°N to 7°N. The 
longitudes covered are 
70° W to 35°W (area W 
or 255), 3° W to 0° 
(area X or 256), 0° to 
35°E (area Y or 257). 


The data are normally 
supplied on a 
latitude-longitude grid 
at 2¥,° resolution. 14,° 
resolution is available if 
required for fields of 
surface pressure, wind 
and precipitation. 


Fields are available at 
6-hour intervals from 
the analysis (T+0) to a 
forecast time of T+36. 
Fields of surface 
pressure, wind and 
precipitation are 
available at 3-hour 
intervals on the finer 
grid. 


5. Back-up arrangements 


Although modern computing systems are very reliable and have operational performance figures that 
approach 100% it is still inevitable that occasions will arise when the computing facilities at Bracknell, 
especially the Cyber 205 used for running the numerical models, are non-operational. At such times 
there is a need for reliable back-up arrangements, so that both chart and coded products can be supplied 
to the user without noticeable interruption. A computer outage of less than 12 hours is covered by 
ensuring that forecast data (in required formats) are routinely prepared for a period of 12 hours beyond 
that normally required. By storing the data on the front-end processing system it is then possible to 
supply operational users on the basis of a 12-hour-old forecast run without the users becoming aware 


of it. 


Regional coarse-mesh products 


Surface pressure, surface wind 
(10 m). 


Heights, winds and 
temperatures at 850, 700, 500, 
400, 300, 250, 200, 150, 

100 mb. 


There are 5 areas. Four span 
25°N to 75°N; the longitudes 
are 90°W to 45°W (area Z or 
205), 45°W to 0° (area N or 
206), 0° to 45°E (area O or 
207), 45°E to 90°E (area P or 
208). The fifth area (area I or 
201) covers 75°N to 90°N and 
90° W to 90°E. 


The data are supplied on a 
latitude-longitude grid with a 
2¥,° resolution. For the polar 
area (area 201), the resolution 
is 2¥,° latitude by 10° 
longitude. 


Fields are available at the 
following forecast times: T+0, 
T+6, T+ 12, T+18, T+24, 
T+30, T+36, T+42, T+48, 
T+60, and T+72. In addition, 
surface pressure and 500 mb 
height are available at T+96 
and 7+ 120. 


Global coarse-mesh products 


Surface pressure, surface wind 
(10 m), 500 mb height. 


Maximum wind and 
tropopause. 


Winds and temperatures at 
850, 700, 500, 400, 300, 250, 
200, 150, 100 mb. 


There are 16 areas. In each 
hemisphere there are two 
polar areas covering latitudes 
70° to 80° and 80° to 90° and 
4 mid-latitude areas spanning 
20° to 70° in quadrants. A 
further 4 areas spanning 20°N 
to 20°S cover the tropical belt. 


The mid-latitude areas are at a 
resolution of 2¥,° latitude by 
5° longitude; the tropical belt 
is at 5° latitude by 5° 
longitude. In the polar regions 
there is reduced longitudinal 
resolution. 


Fields are available at 6-hour 
intervals from the analysis 
(T+0) to a forecast time of 
T+48. In addition surface 
pressure and 500 mb height 
are available at 7+60, 7+72, 
T+96, and 7+120. 
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To ensure coverage beyond 12 hours there is a mutual exchange of global forecast data between 
Bracknell and Washington, using GRIB code as a medium. Data from Washington arrive at Bracknell 
too late to be used operationally in real time, but after decoding and storage in a compatible output data 
set they can be used to give an almost complete range of products, although again from a 12-hour-old 
forecast run. 

The obvious omissions from the product list when operating in a back-up mode are those from the 
fine-mesh NWP model. Although in the recent past it has proved possible to run versions of the previous 
Bracknell NWP model on the front-end processors, this is now considered to be increasingly 
unsatisfactory and other sources of high-resolution forecast data are required. It is possible that other 
members of Revion VI will be able to supply such data for back-up purposes. 
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